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Abstract 



The objective of reliability sensitivity analysis is to determine input variables that mostly contribute to the variability 

of the failure probability. In this paper, we study a recently introduced method for the reliability sensitivity analysis 

based on a perturbation of the original probability distribution of the input variables. The objective is to determine 

the most influential input variables and to analyze their impact on the failure probability. We propose a moment 

independent sensitivity measure that is based on a perturbation of the original probability density independently for 

t— I each input variable. The variables providing the highest variation of the original failure probability are settled to be 

more influential. These variables will need a proper characterization in terms of uncertainty. The method is intended 

r_^ to work in applications involving a computationally expensive simulation code for evaluating the failure probability 

r*\ such as the CO2 storage risk analysis. An application of the method to a synthetic CO2 storage case study is provided 

together with some analytical examples. 

Keywords: sensitivity analysis, reliability analysis, uncertainty analysis, failure probability 



1. Introduction 

Carbon Capture and Storage (CCS) stands for the collection of CO2 from industrial sources and its injection 



> 

^j into deep geological formations for a permanent storage. There are three possible sites for injection: unmined coalbed 

formations, saline aquifers and depleted oil and gas reservoirs JT]. Nevertheless, the following principal environmental 
question arises: what is the probability that CO2 will remain underground for hundreds to thousands of years after its 

—J capture and injection into a storage formation? 

The primary risk of CO2 geological storage is unintended gas leakage from the storage reservoir J2|[3). In this 
work, we focus on a leakage from the storage formation through a fault or a fracture. This can happen when the 
reservoir pressure is higher than the caprock fracture pressure. Numerical modelling and simulation has become an 
J> integral component of CO2 storage assessment and monitoring. The reservoir simulation models are constructed based 

on the reservoir and production data. They are used to predict and analyze CO2 plume distribution and the reservoir 
pressure development during the injection and the storage periods. With the help of the numerical simulation models 
it is possible to forecast CO2 storage performance and to evaluate the risks of a possible leakage. 

Risk and uncertainty analysis has been recognized as a principal part of safety and risk assessment. Generally 
speaking, when the main sources of uncertainty have been identified, Uncertainty Analysis (UA) is focused on quanti- 
fying the uncertainty in the model output resulting from uncertainty in the model inputs. At the same time, Sensitivity 
Analysis (SA) aims to identify the contributions of each model uncertain input to the variability of the model out- 
put. Uncertainty analysis may be equally performed to assess the reliability of the system. A typical example of a 
failure probability estimation in the CO2 storage risk analysis is the estimation of the probability of exceeding the 
caprock fracturing pressure during the CO2 injection phase. If we denote P reS ervoir the reservoir pressure and P 'fracture 
the caprock fracturing pressure, then we consider the following failure probability: 

Pf — ryr reservoir — ^ fracture)- 
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In practice, uncertainty and sensitivity analysis require a large number of reservoir simulator runs to explore all the 
input variables space. However, higher accuracy of a simulator usually results in a higher simulation time. One 
simulator run can take from few minutes up to several hours or even days. Therefore, when the simulation time 
becomes too high, uncertainty analysis may become unfeasible. For this reason, in this work, in order to estimate the 
failure probability pf for the expensive reservoir simulator, we use Gaussian Process (GP) response surface model, 
also known as kriging |@1[5]|6). The kriging method was originally introduced in the field of geostatistics by Krige 
in the 1950's |7| and formalized in 1960's by Matheron (8). In (4) Sacks et al. proposed the statistical approach to 
uncertainty analysis of complex computer codes referred to as the Design of Computer Experiments. In a nutshell, the 
approach consists in building an approximation of the reservoir simulator input/output relationship starting from a set 
of simulation runs at a carefully chosen input variables configurations referred to as experimental design or training 
set. The obtained response surface model can then be used to predict the model output for a new non simulated input 
with a negligible computational time. Therefore, uncertainty and sensitivity analysis become affordable. 

One of the most challenging problems in risk analysis is to identify the failure region and to compute the failure 
probability. However, the failure probability will usually depend strongly on the probability distribution of the input 
variables. Reliability Sensitivity Analysis can help in understanding the relationship between each input variable 
uncertainty and the failure probability. The problem is to identify the set of input variables that need to be well 
characterized in terms of uncertainty distribution. 

For the time being, numerous approximation and simulation methods are available for estimating the failure prob- 
ability such as First / Second Order Reliability Methods, Monte Carlo sampling, importance sampling, directional 
sampling, subset simulation, etc. (9j [TUJ Qj] [12] . However, there are few sensitivity analysis methods developed for 
the failure probability analysis. The widely used methods for sensitivity analysis are based on a variance decomposi- 
tion of the output. Given the probability distribution of the input variables, Sobol indices are expressed by the ratio of 
the variance due to a given input on the total output variance fOl . Knowing the probability distribution of the input 
variables and the output, we can so define sensitivity indices for each of the input variables. Variance based indices 
are usually interesting for measuring the input/output sensitivity, however they can be poorly relevant to our problem 
of evaluating the impact on failure probability. 

There have been few attempts to develop sensitivity analysis methods well suited for reliability analysis. First, 
as complementary results of the First Order Reliability Method, sensitivity to the distribution of the input variables 
can be obtained. Sensitivity is expressed as the partial derivative of the reliability index /3 |9|. Another approach was 
proposed by Morio 1 14]. Therein, the author uses the variance decomposition and Sobol' sensitivity indices to study 
the rate of change in the failure probability due to the changes of the input distribution density parameters. Borgonovo 
et al. [15] suggested some moment independent importance measure in the reliability analysis. This measure does not 
involve the variance. For a fixed variable jc, it quantifies the effect of knowing x, by computing the L\ norm between 
the unconditional joint density f s (-) and the conditional density /j|. v ,(-)- 

In this paper, we study a moment independent approach for sensitivity analysis of a failure probability |[T6l . 
The influence of the input variables on the failure probability is obtained by perturbing the prior probability density 
function fi{x). In particular, we estimate the effect of the perturbation on the value of the failure probability p/. 
Here, we propose to distinguish distributions classes by their supports. For the case of a bounded support, such as the 
uniform or the triangular distributions, the main source of uncertainty is about the boundaries of the support. On the 
other hand, in the case of infinite support, such as normal or log-normal distributions, the main source of uncertainty 
comes from the distribution parameters, such as mean and variance. The estimation method has the advantage of 
being very efficient in terms of number of simulator calls. In order to estimate the sensitivity indices for all the input 
variables, the performance function is evaluated only once on a Monte Carlo sample used to estimate the reference 
failure probability pj. 

Our paper is organized as follows. First, we introduce some density perturbations for different families of dis- 
tributions. Later, we introduce the technique to compute a perturbed failure probability using the same Monte Carlo 
sample. This is based on an inverse importance sampling technique [ 17 1. Finally, we present the formulation for the 
moment independent sensitivity indices and demonstrate its applicability on an analytical and a CO2 storage reservoir 
case examples. 



2. Density perturbation influence to failure probability 

Let us denote by g(x) the performance function of the system, x e Q c Mr is a set of independent input variables 
with the joint density f x (x) - Yltsi fxfod- The failure probability is expressed as: 

Pf = P(g(x) < 0) = E A [l^)< l = f /*(J0<tt, 



where Qy • = {x e Q : g(x) < 0} is the failure region. 

In general, the distribution density f x (-) is provided by experts on the basis of some indirect measurements or 
some limited observation data f!8l . Here, we study how a perturbation of the original probability density /*(•) affects 
the failure probability of the system pf. We assume that the input variables x,, / = 1, . . . ,d are independent random 
variables with marginal densities f Xr So that, f x (x) - n,=i fx,(xi). 

This work is inspired by previous work |[T8l[T6l . Originally the method was applicable only to the normal distri- 
butions. Here, the method is extended to the more general case of an exponential family. The objective of this study 
is to estimate the influence of the input random variable from the vector x to the failure probability of the system pj. 
The proposed approach consists in perturbing the original density for a given fixed variable Xi while keeping constant 
the probability density functions for all the other variables x-t = (xi, . . . , Xi-\, Xi+\, . . . , Xa). Then, a new value for the 
failure probability is computed. If this new value p^ differs significantly from the reference value pj, it means that 
this selected input variable x, is influential. Conversely, if the new failure probability p is is close to pr, then the input 
xi has low influence on the failure probability. 

2.1. Density perturbation for an exponential family 

To define the probability density perturbation, first recall the definition of an exponential family. An exponential 
family is a set of distribution having density function that can be expressed in the form of: 

f x {x\0) = h(x) exp (iKQTGt) - A(0)) , (1) 

where is a vector of distribution parameters, tj(0) is a natural parameter, T(x) is a vector of sufficient statistics, h(x) 
is an underlying weight function and A(0) is the cumulant generating function [19]. The cumulant generating function 
ensures that the distribution integrates to one, i.e.: 



/ 



A(0) = log h{x)exp(r](0)T(x))p(dx), 

where p. is the reference measure (for example Lebesgue measure). It could be continuous or discrete. In this paper, 
we mainly consider continuous measure. However, all the calculations are valid for a discrete measure as well. Expo- 
nential family contains most of the standard discrete and continuous distributions that we use for practical modelling, 
such as the normal, Poisson, Binomial, exponential, Gamma, multivariate normal, etc. 

To define the density perturbation for this family of distributions, we use some ideas coming from information 
theory ll20ll . Kullback-Leibler (KL) divergence is used to measure the magnitude of a perturbation. KL divergence 
quantifies the "closeness" of two probability distribution P and Q. Suppose that P and Q are continuous probability 
distributions with densities p(x) and q{x) (with respect to Lebesgue measure). Then, the KL divergence between P 
and Q is given by: 

D KL(^0= f Pimog^rdx (2) 

J-™ q(x) 

For 6 > 0, let us denote for the variable x, the new (perturbed) density as f Xir {-),- We select f XlT (-) in such way that: 

DKL(/*./*) = & (3) 

Possible values of the perturbation 6 may be restricted by some inequalities on Kullback-Leibler divergence ED . If 
we define the function r(x) := t§,(* e Q c R d ) and assume that < r < r(x) < R for all x e £2 c R d . Then, 



according to 12T1 . we have: 

Dki>.<?) * ° 

DkL(P.9) ^ 4rjg = <W 

According to these inequalities, we choose 6 e [0,S majc ], where 5 max - 4 ~^ can be computed precisely. 

Let us consider the original density f Xi (x) = h(x) exp{rj(6)T{x) - A(ff)) from exponential family. In order to stay in 
the same family, we propose to restrict the choice of possible perturbations among the following class of densities: 

f XiT (x) = exp(rT(x) - ^{r))f Xi {x) = h(x) exp (T(x)(t}(0) - t) - (A(ff) + ^t))) , (4) 

Here, t is a constant depending on 5 (it is chosen under the condition <j3j) . The function ij/if) is a normalization 
function and it may be expressed as: 



\f/{T) = log 



f 

•J — Cx 



exp( T T( t))f Xi (t)dt 



It is the cumulant generating function for the perturbed probability distribution f Xir (x), Moreover, if yU and a 2 are the 
mean and the variance of the original probability distribution, then: 

«A(0) = 

<A'(0) = ix (5) 

<T(0) = O- 2 

We aim to perturb f Xi in such a way that the KL divergence between the original density /„■ and the perturbed density 
f XjT is equal to 6. Notice that: 

X°° f (t) r°° 

f XiT (t)\og J -f^-dt = f XiT (t)(rl(t)-^(T))dt (6) 

= TiA'(r) - «A(t). 

Hence, t should satisfy the equation: 

Tlf/if) - iff{T) = 6. (7) 

Let t* = t((5) be one solution of (|7J. We use this parameter in order to define the perturbed density modification f XiT (-) 
defined by Q. 

Now, let us consider the function: 

G(T) = Tlf,'(T) - <KT) " 5. 

This function has a global minimum at r — 0: G'(t)| t=0 = 0, G"(t)| t=0 = iff"(r) > and G(0) = -5 < for d > 0. 
Moreover, G'(t) = ti//"(t): G'(t) < 0, (t < 0) and G'(t) > 0, (t > 0). Thus, G(t) is strictly decreasing f or t < and 
G(t) is strictly increasing for r > 0. Hence, the function G(t) has not more than two zeros T\ < and T2 > 0, if both 
of them 7"i and T2 fall into domain of the function i/Kt). 

For every fixed level of 6, we can study two possible effects of the perturbation Q. We denote the corresponding 
perturbed densities by f Xjr and f Xi ^ . Then, the joint perturbed probability density is expressed as: 



fx^m = f xlT . f] f xt (x k ), 7=1,2. 



yfc=l,jt#i 

The corresponding value of the perturbed failure probability p is (j — 1 , 2) can be computed as the following integral: 



P®j = %%. [*«(*)«)] = I ^xo/ivr/^' ./' = 1,2. 



(8) 



In the same way, the interaction effect can be estimated by perturbing two variables xt and Xj at the same time by 5\ 
and 62 respectively. Suppose that t(6\) = (ti^.t^,) andr^) = (T\g„T2s 2 ) are the solutions of equation (|7), where 61 
and 62 are the perturbations of KL divergence (|3]l for the variables xt and Xj, respectively. The new joint probability 
density function is: 

d 

Jxij,AS\).T{s 2 y X > ~ J x M.i\)J x Mh.) I I fxk( x k)- (") 

k=l,k*i,j 

The corresponding value of the perturbed failure probability is estimated in the same way by putting in ([8]) the new 
joint probability density / %rHl)T( , 2) (0- 

In the next section, we introduce a method to estimate efficiently the perturbed failure probability pa using the 
same Monte Carlo sample as for the estimation of the original failure probability pf. First, we study the effect of the 
perturbation for different probability distributions. 

2.2. Resulting distributions 

Here, we provide a summary table for the considered distributions from exponential family with the resulting 
perturbed distribution. For the case when the natural parameter rj(6) is a vector of functions (like normal and log- 
normal distributions), we propose to analyze the perturbation effect separately for each of the components of J](6). For 
example, for the normal distribution: 



T(x) 



x 

.-2 



T](n,cr) 



p/0- 2 

-l/2cr 2 



Then we propose to analyze two different density perturbations Q by taking: 



Table [T] provides the results obtained for some distributions from the exponential family: normal, log-normal, expo- 
nential and Poisson (as discussed all the calculations are valid for a discrete probability measure as well). As it can 
be seen, the new perturbed density f XjT is still in the same family of distributions with the perturbed distribution pa- 
rameters. For the case of normal and log-normal distributions when we study the effect of perturbation of the second 
component of the natural parameter, there is only one solution for t: 1 - 2rcr 2 > 0. This solution could not be found 
analytically but with the help of a numerical solver. For the exponential and the Poisson distributions the solutions 
for t is expressed with the Lambert W function. This stands for the the multivalued inverse relation of the function 
f(w) - wexp(w), where w is complex. We denote by Wq(x) the upper real branch of the Lambert function on the 
interval [- 1 /e, 0] and by W- 1 (x) the lower real branch on the same interval. 

We will analyze the effect of these perturbations on an analytical example in Section[6] 



Distribution 



Natural Parameter 
10) 



Sufficient Statistics 
T(x) 




New cumulant function 
Mr) 



t\ (6) and 
t 2 (<5) 



Resulting distribution 



Normal NQj,cr 2 ) 

LogNormal iogNip, a 2 ) 

Exponential Exp(/1) 
Poisson Pois(/t) 



-l/2o- 2 

H/0- 2 
-l/2o- 2 

-A. 

log(i) 



A«"+ — 
^ + 1 log(l - 2to- 2 ) 

j^, + i log(l - 2tO 
logfe) 

<«e r -l) 



Numerical solver 



Numerical solver 



Tl,2(5) = 



A(w^„(-e-'-')+i) 



Tl,2(<5) = W-l,o(-^) + l 



r Nip + to- 2 , 0- 2 ) 1 

logAA(/j + to 2 , a 2 ) 

Exp(,l - t) 
Pois(/i exp(r)) 



Table 1: Resulting distributions. 



3. Bounded support distribution 

Now, let us consider xt ~ U[a, b], the uniform distribution on the interval [a, b]. The density is expressed as: 

fxi(x) = : he[a,b](x),(b > a). 

b - a 

The uniform distribution does not belong to the family of exponential distributions. It has a limited and a bounded 

support [a, b]. If we apply the same perturbation as for exponential families, the normalization function becomes: 

(10) 
\ t\d - a) j 

Then the equation for t(<5) is: 





1 Tb _ ™\ 




W\ T > IU 6 // , ,TLM. 




\ lib - a) ) 


Tbe Tb 


- Tae Ta - e Th + e™ / e Tb - e™ 




e h - e™ * 6 1 r(& - a) 



This equation has no explicit solutions for r. The solutions can be found using a numeric solver. Suppose that 
t* = t((5) is a solution of equation ( fT0] i. Then, the perturbed density is: 

T*e T *- v 
Therefore, the new perturbed variable jc, t is no longer uniform on [a, 6]. This density modification for a = — 1, Z? = 1 



and 6 - 0.5 is displayed in Figure 1(a) 



Notice that working with uncertain variables defined on a compact support, the main source of uncertainty is on 
the boundaries of the support. For such distributions with a bounded support, we propose to apply another density 
perturbation. The idea consists in perturbing the original boundaries by r = ±<5. In the same way as with infinite 
support we consider the effect of positive or negative perturbation. For example, consider x-, ~ U[a, b] to be uni- 
formly distributed on the interval [a,b]. Then, in order to stay inside the support the perturbed random variable x iT 
is uniformly distributed either on U[a + 6,b] or on U[a,b - 6]. The corresponding density for perturbed uniform 
distribution can be expressed as: 

1 1 

frl(x) = -he[a+d,b]( x ) or M(x) = zIxE[a,b-S](X). 

b — a — o b — a — o 

The same perturbation may be applied to a triangular or a trapezoidal distribution. It can be also applied to the 
truncated Gaussian distribution if one is interested about the boundary influence on the failure probability pf. In this 
case the density function should be corrected for the new boundaries. 



Uniform density perturbation 



Uniform density perturbation 










! ! 


























Original density f(x) 
















Negative t 
Positive z 

















(a) Exponential density modification 



(b) Boundaries perturbation 



Figure 1 : Uniform density perturbation. 

Next, we explain how to estimate efficiently a perturbed failure probability pf 6 with no additional CPU cost. 



4. Inverse importance sampling and sensitivity analysis 



Monte Carlo sampling is one of the most popular simulation methods to estimate a failure probability. We consider 

''. Recall that all the input variables are independent and that f x (x) = Y\k=i fx k ( x k) ls 

ff , _ _ , i.i.d. 



the input variables space Q e 

the joint density of the input variables. Let X 
the failure probability pf is given by: 



1*1. 



< X N1 



/j(-) be a sample of size N. Then, the estimation of 



Pf 



1 N 



gm<o- 



(ID 



Now, assume that f Xjr is a new perturbed density for the input variable x,-. Then the new joint density is f Xjs = 
fxi ■ ■ ■ fxi-Jx-JxM ■ ■ ■ fx d = fx, T (xd n^i,^, fx k (Xk). The corresponding failure probability p i6 is defined as an expecta- 
tion of the indicator function: 






Here, we propose to apply the technique used in the Importance Sampling (IS) simulation method. We multiply the 
integrand function by 1 = 444| . Both density functions fx(x) and fj- is are the products of the density functions of the 
independent variables x - (xi, . . . ,Xd) eficl 1 ' with the only difference for the variable x,. Therefore, we obtain: 



r r fx oa) 

Pis - l g (x)<ofx k ,(x)dx - I g (x)<o " , Mx)dx = E f . 
Jn Jn Jxi(Xi) 



g(x)<0 



fx t (Xi) 



By doing so, we do not need to throw a new sample according to the unknown density function fx iS {x). We are 
working in the same probability space integrating the function 



I, 



f HT m 



probability pa we keep the same sample points from the failure region: X'Y = \x e X^ : g(x) < 0} that provide 



. So that, to estimate the perturbed failure 

N . 



f 

non-zero values of the indicator function I g (i)<o- The estimation of the failure probability for the perturbed density is 
expressed as: 



k=\ 



)<() 



fx, T (Xki) 
fxi(Xki) ' 



(12) 



If we are interested in the interaction effects, we perturb the probability densities for the variables xt and Xj simulta- 
neously. Then, the new joint density is: 

d 

k=\Mi,j 

Therefore, in this case the new failure probability can be estimated by: 

— 1 y T Jx iT( s l} \ x ki)jXj 



(%) 

This describes the interaction effects of two variables x, and Xj on the failure probability pf. 

The proposed reliability sensitivity analysis is based on the analysis of the value of the perturbed failure probability. 
In order to clearly differentiate the magnitude of the influence to the failure probability, we propose a sensitivity 
measure for the input variables. We provide the indices formulation in the next section and we study their statistical 
properties in | Appendix B| 

5. Sensitivity indices formulation 

There are many possible choices of sensitivity indices. In this paper, we propose the sensitivity indices based on 
the difference p^ - pj with the original failure probability pj. It is expressed as ratio: 

s*-»^*L (13) 

Pf 

Another possible formulation could be found in Lemaitre and Sergienko et. al, 2012 [ 16]. 

The support of Sts formulated in ( p"3] > is [-l,+oo). A negative value of Sis means that the proposed density 
modification reduces the failure probability. Conversely, a positive value of this index means an increase in the failure 
probability. Zero value of S a means that the variable Xi has no impact on the failure probability. 

In practice, pf is estimated by p~f with the Monte Carlo simulation method. In the same way, according to ([12} p is 
is estimated by ~p j6 . The estimator of the indices S ,s can be expressed as: 

^ Pis~ Pf 

J iS - — • 

Pf 

This estimation provides an asymptotic unbiased estimation of S is- Moreover, according to the Central Limit Theorem 
(CLT) and the A - method we have: 

(Su-Sis) -» JV«U) 



V var aN 



The proof and an expression of the asymptotic variance can be found in |Appendix B.2| 

The above basic indices are computed separately for every fixed value of the density perturbation 5. It allows to 
study the effect of perturbation by varying the value of 6. It can help in reliability design optimization by adjusting 
the distribution parameters of the input variables x = (*i, . . . , Xj) in order to achieve the lowest failure probability pf, 

6. Analytical function example 

To investigate the previous indices, we consider a linear function of three independent normally distributed vari- 
ables: Xi,X2,X3 ~ N(Q, 1): 

g(x u x 2 ,x 3 ) - 3 -O.Lti -0.5x2 - l-0x 3 . (14) 

S 



As a linear combination of independent Gaussian random variables, g(x) is distributed as N(fi, &), where // = 3 and 
<x 2 — 0.1 2 + 0.5 2 + 1.0 2 = 1.26. The original failure probability can be explicitly computed: 



p f = 1 - d>(^-) = 3.7 x 1(T 3 . 



The estimation provided by a Monte Carlo sample of size N = 10 6 yields p~f = 3.69 X 10 3 . We use the same Monte 
Carlo sample to estimate the perturbed failure probabilities p is and the sensitivity indices Sis. 
Recall, that the natural parameter of the normal distribution is given by the vector: 



770", o") = 



-l/2cr 2 



First, we consider the effect of the perturbation of the first component of natural parameter 77. This perturbation refers 
as a shifting of the original mean of the distribution (see Table fT). Figure Pf] depicts sensitivity indices calculated for 
6 € [0, 1] with the positive ( |2(b)[ ) and the negative ( |2(a)| i values of perturbation parameter r. It can be clearly observed 
that the highest impact on the failure probability is due to the variable X3 and that the variable x\ has the lowest impact 
for both cases of t. Moreover, the higher the value of S, the higher the influence of the variables. 
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(a) Negative perturbation 



(b) Positive perturbation 



Figure 2: Basic sensitivity indices example. Mean shifting. 

Now, we study the effect of the perturbation of the second component of the vector of natural parameters -l/2cr 2 . 
In this case, we perturb both the mean and the variance of the original distribution (see TablefTJ). We use a numerical 
solver to find the solutions for r. In this case we consider the only possible effect of perturbation when the normal- 
ization function ip{j) is defined and 1 - 2tct 2 > 0. Figurepldisplays the estimated sensitivity indices in this case. We 
can observe the same relationship: the variable X3 has the highest influence on the failure probability. However, the 
magnitude of the influence differs from the case where we only consider a mean shifting. The variable x\ still has 
almost negligible impact. 
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Figure 3: Basic sensitivity indices example. Mean and variance shifting. 

7. C02 storage case example 

In this section, we consider a CO2 storage reservoir simulation model. As discussed in the introduction, subsurface 
CO2 storage is always associated with an excess reservoir pressure. One of the primary environmental risks is a 
pressure-driven leakage of CO2 from the storage formation. 

In order to assess the risk of CO2 leakage through the cap rock we consider a synthetic reservoir model. The 
structure of the reservoir is reduced to its simplest expression. The model is made up of three zones (Figure Q): a 
reservoir made of 10 layers, a cap-rock made up of 1 layer, a zone-to-surface composed of 9 layers. 




Figure 4: Reservoir model. 

The XY size of the grid is set at 10 km total length representing 26x26x20 model grid. Each layer is 5m thick, 
including the cell above the cap-rock. The zone above the cap-rock (up to the surface) is currently set to 1 layer. The 
salinity of the water is 35gm/l. The temperature of the reservoir is set to 60C and the initial pressure is hydrostatic. 
The injection bottom rate is set to 10 6 tons/year. The fracture pressure is estimated by geomechanical experts to be 
P 'fracture = 122 bars. Exceeding this value during the injection can lead to a leakage. The simulation period is 55 years 
that include an injection period of 15 years followed by 40 years of storage. In this study we analyze the possibility 
of a leakage through a cap rock. Therefore, we consider pressure in the storage reservoir at the last year of injection 
as an objective function. 

The uncertain variables selected for this study characterize the reservoir and the fluid properties. It implies different 
CO2 flowing possibilities between the reservoir layers. Table (|2} represents the variables description with their range 
of minimum and maximum values. 
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Name 



Description 



Min Max 



PORO 

KSAND 

KRSAND 



Reservoir Porosity 0.15 0.35 

Reservoir Permeability 10 300 

Water relative permeability end-point 0.5 1.0 



Table 2: Uncertain variables. 



For sake of clarity we transform the original intervals into [—1,1]. In this section, we assume the truncated standard 
normal distribution for all the input variables A/[-.i,i](0, 1). 

The performance function for this example can be formulated. Suppose that V reservo ir(x) is a function of the 
reservoir pressure depending on the input variables configuration x e Q c R d . Then, the performance function 
defining the event of a gas leakage can be expressed as: 



g(x) = Pfra. 



-P. 



reservoir 



(X). 



The reservoir pressure P rereno ,>(x) is computed with a complex dynamic reservoir simulator. For this reason, we use 
the Gaussian process based response surface model approximation P(x) QH]|5]|6|. By approximating the function of 
the reservoir pressure, we can quantify the risk and estimate the reliability of the system. 

The original reference value of the failure probability computed by GP model approximation with Monte Carlo 
sample of size N = 10 6 provides an estimation: p~f - 2.26 x 10~ 4 . We keep this sample to estimate the perturbed 
failure probability p^ and the sensitivity indices Sjg. 

There are two possible ways to perturb the truncated Gaussian distribution. We can use the perturbation defined 
for the exponential family or we can study the effect of the perturbation of the distribution support boundaries. In this 
example, we compare the results for both cases. We start with the sensitivity indices calculated by the mean shifting. 
Figure |5]displays the evolution of the sensitivity indices for 6 e [0, 1] for negative (Figure 5(a) i and positive (Figure 
5(b)| > values of shifting. The variables ranking depends on the sign of r. When r > (i.e. the positive mean shifting) 
the porosity variable has the highest impact on the failure probability. It means that increasing the mean value of the 
reservoir porosity PORO leads to increasing the failure probability. On the contrary, increasing the mean value of the 
reservoir permeability KSAND and the water relative permeability end-point KRSAND has a negative effect on the 
risk of leakage. For the negative mean shifting, the variables KSAND and KRSAND have the highest influence on 
the failure probability. Reducing the reservoir permeability and the end-point water relative permeability impedes the 
gas flow in the reservoir. It increases the risk of the leakage. 



Sensitivity Indices 



Sensitivity Indices 
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Figure 5: Mean shifting. 
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Now, we consider the boundaries perturbation. It means that we are moving one of the the distribution boundaries 
by 6 in positive or negative directions keeping the values of mean and variance unchanged. Figure [6] depicts the 
sensitivity indices for 6 e [0, 1] for negative (Figure [6(a)| and positive (Figure [6(b)] > values of t — ±5. 
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Figure 6: Boundaries shifting. 

When t < (i.e. the resulting distribution is N\-\,\-gtf$S, 1)) the porosity variable has the highest impact on the failure 
probability. It means that by decreasing the maximum value of the reservoir porosity PORO the failure probability 

It can be also observed that for this variable the new perturbed failure 
When r > the resulting distribution is A/[-i+<s,i](0, 1). For this case, 



decreases. It is also shown by Figure 5(a) 

probability pig is equal to zero when 5 > 0.3. 

increasing the reservoir permeability KSAND and the water relative permeability end-point KRSAND reduces the 

failure probability ptg. 

Both methods provide comprehensive and complementary results. If the main uncertainty is about the boundaries 
the one can start with the boundaries perturbation. By moving the boundaries of the original distribution, the one can 
determine the safe intervals for the input variables by detecting the value of S: pi$ = or S , # - -1. After that, the 
effect of the mean perturbation can be studied. 



8. Conclusions 

In this paper, we have studied and adapted a recently introduced approach to the reliability sensitivity analysis. 
Currently the majority of the methods for reliability analysis is based on the variance decomposition and SoboF 
sensitivity indices. We present a moment independent sensitivity measure. The method is based on a perturbation of 
the original probability distribution of the input random variables. We can analyze the a priori assumption about the 
input distributions and measure the effect of some possible deviations from this assumption. In particular, we select 
the Kullback-Leibler divergence as a measure of the perturbation. 

We have provided different possible density perturbations with the resulting distributions for the exponential fam- 
ily of distributions. We have also studied the distributions with a bounded support and the effect of the boundaries 
perturbation. 

Considering a proposed perturbation for an input variable, we present an effective method to estimate the corre- 
sponding perturbed failure probability. The method is based on a technique coming from the importance sampling 
simulation method. It allows to estimate the new failure probability without supplementary performance function 
evaluations. The new sensitivity indices formulation describes the relationship between the new failure probability p/g 
and the original failure probability of the system pj. By varying the value of the perturbation 6, we can study how the 
positive or negative probability density perturbation affects the failure probability. If the model has controllable input 
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variables, the method can help improving the system reliability and the design optimization. The presented analysis on 
the statistical properties of the proposed estimators for the perturbed failure probability ~pa and the sensitivity indices 
S is shows asymptotic normality of the estimators. 

We investigated the method on an analytical and a CO2 storage reservoir cases. The method provides promising 
results and can be applied in the reliability sensitivity analysis. 
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Appendix A. Properties of the normalization function 

For a given random variable x with probability density f(x). We propose the density modification x T ~ / T (-) as 

follows: 

f T (x) = exp(rx - (A(t))/(x), 

where if/(r) is a normalization function given by: 



iff(r) = log 



f 



exp(Tx)f(x)dx 



Let us define T> — {t e R : i^(t) < +00}. £) defines the interior of T). We will also suppose, that 3e : T> z>] - e, s[. 
Here, we will study the properties of this normalization function. 



lendix A.l. Derivatives 




• f(T) = E[x r ],TE0 


Tt [£L exp(Tx)/(x)rfxl r°° 
W(T) = roo = 1 xexp(rx 
J exp(Tx)f(x)dx J-00 


• i//'(t) = E [x T - E(x T )] 2 , t e D 


d 

<A"(r) = - 
dr 


/-*oo 

I x exp(rx - if/(T))f(x)dx 

<J — oo 



Xfr 

<J — CO 



(x)dx = E[jc t ] 



[x(x - if/'(r)) exp(rx - \p(r))f(x)dx\ - 

x f T (x)dx - i//(t) I xf T (x)dx = I x f T (x)dx - [E[x T ]] = 

CO v/— CO %J — CO 



= E[x T - E(x T )] 2 = VAR(x T ) 

n 3 



T (x)dx ■■ 



if/'"{T) = E [x T - E(x T )Y ,reD 

Xco 
[x->I/(t)TM 
CO 

XCO /-»CO y^»CO 

x 3 f T (x)dx - 3i//(t) I x 2 f T (x)dx + 3iA'(t) 2 I xf T (x)dx - i/r'(r) 3 
00 <J — oo %J —CO 

XCO /-»CO 

x 3 f T (x)dx - 3iff'(r) I x 2 f T (x)dx + 2i//'(t) 3 = if/"(r) 
00 *J — oo 
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Appendix B. Statistical properties of the indices estimator 

Here, we will study statistical properties of the estimator of the perturbed failure probability ~p~js and corresponding 
estimator of the sensitivity indices S is - P '~ — ■ We will start with studying the properties of /T^. 

Appendix B.l. Estimator of the perturbed failure probability 

Suppose, /x(x) = nf=i fxi(xd is the input joint density and f XjT is a perturbed probability density for the variable 
X[. Recall that for a sample of size N: {Xi, . . . , X#} ~ ' / x (x), the estimation of ~p is is computed by: 

(V 



- _ 1 Vt /v, r (Xt,) 

Pit ~ n 2-t 1 i' (x ' < °)" 



*£-} '^ /,(*>) 
First, we study the expectation and the variance of this estimator. 

1. E fl \pis\ =p iS 

2. VAR /x [p is ] = lVAR /s 



T /v,y(-r<> 



This variance tends to when N — > oo. Furthermore, by the Central Limit Theorem (CLT): 

1 



J IjKxxo J^/x(x)dx - p: 
)y the Central Limit r 

(Pis - P«) ^ AT(0, 1) 



VVAR/ X [£, 

Note that the covariance between the estimator p~f and p^ does not vanish. Indeed, we use the same sample to 
estimate pf and p is . We can compute this covariance: 

COW (pf,p iS ) = E f% (p/p iS )-E(p f )E(p is ) = 

= jjPaO- - Pf) 

The value of this covariance decreases when the sample size N increases. 

Appendix B.2. Sensitivity indices estimator 
Recall, that the first sensitivity index is: 

„ Pis ~ Pf pis . . . , 

Sig = = l,i = 1,.. . ,d. 

Pf Pf 

We estimate this value with the Monte Carlo method and the importance sampling by estimating consistently pig and 

Pf. The estimator of this index is: _ 

-jr Pis . 

Pf 

Here, we will study some proprieties of this estimator. 

It is not straightforward to compute directly E/ s \S i$\ - E/ x W- —\ and VARy x \S -A = VAR/ s |£ . We propose 

to use the Delta Method to approximate these values ll22ll . 

Let us recall the Taylor expansion with integral form for the remainder. Let <p be a two times differentiable function 
on [to,t], then: 

(f>(t):=<t>(t ) + <p'(t )(t-t )+ ( (1 - u)<p"(u)du (B.l) 

J to 



We will define a function: 



MA y(t) 1 

4>(t) = -T- - 1, 

x(t) 



where x(t) — (1 - f)p/ + tpf and y(0 = (1 — f)p« + tpis- For this function: 0(0) = S,<5 and 0(1) = Sis. Following the 
Taylor expansion ( |B.1[ ) we will expand 0(f) with f = 1 and f ( > = 0. First, we will compute the derivatives. 
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1. x'(t) - Pf- Pf 






2. y'(t) = J) is - pit 






Then, 


cf>'(t) -- 


y'(t)x(t) - x'(t)y(t) 
xHt) 




^'( f )L -- 


PfPiS ~ PfPiS 



((l-t)p f + tp f f 



The second derivative is: 

2x(t)x'(t)(p f p is - p/Pis) 2x'(t){j}fp i6 - p/pis) 2(p f - Pf)(pfPis- PfPis) 
^ ~ xHf) ~ ~xHt) ~ xHt) ' 

Therefore, the reminder is: 



X 



i 2(1 - 00/ - Pf) [p/Pis - PfPis) _(p f - PfKpfPts - PfPts) 



3 dt - „2, 



((l-t) Pf + tp f ) PfPf 

So that, by the Taylor expansion ( |B.1| > we obtain: 

PfPiS ~ PfPiS (Pf ~ Pf)(PfPiS ~ PfPis) 



1 1<5 — J iS ' 



Pf PfPf 



The last term R - — — p ,_" > PlP '" is the remainder. This remainder is bounded and we can neglect it for the 

r f Pf b 

approximation. — — 

■^ c ^ PfPis ~ PfPiS 
OiS ~ o iS H S • 

Pf 

Now, we can approximate the mean and the variance of S $. 

1. E A \Su]~S u 

2. VAR /x [S is ] ~ VAR /x [|] + VAR /x [^f] - ^COVfe,?/) = ^VAR /x [p is ] - ^^. 

Therefore, the variance of the indices estimator tends to when Af — > oo. With some extra computations we can show 

that: 1 /- v^oo 

(Sis-Sit) -» 7V(0,1). 
Knowing the variance of the estimator, the confidence region for the indices may be computed. 
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